/*****************************************************************************
  File Name      : "phi-GRAPE-ser.c"             // Stellar Cluster version...
  Contents       : N-body code with integration by individual block time step 
                 : together with the using of GRAPE6/GRAPE6a board.
                 :
                 : Normalization to the physical units!!!
                 :
                 : External Potential added 
                 : Plummer-Kuzmin: Bulge, Disk, Halo
	         :
                 : Rebuced to the Single BH -> Plummer !!!
                 :
                 : Stellar evolution added
                 : Stellar lifetimes: Raiteri, Villata & Navarro (1996)
                 : IMS mass loss: van den Hoeg & Groenewegen (1997)
                 :
                 : STARDISK: Drag force... Chingis, Denis & Max... 
	         :		 
                 : STARDESTR: Tidal disruption of stars by BH...  
                 :
                 : Define the SC DC for "a" & "adot" better accuracy
                 :		 
  Coded by       : Peter Berczik
  Version number : 2.1.PUB
  Last redaction : 22.09.2008 19:17
*****************************************************************************/

//#define DEBUG

#define EXTPOT		// external BH 

#define STARDISK	// gas disk around BH + drag
#define STARDESTR	// disruption of stars by BH tidal forces

//#define NORM		// possible physical normalization...
//#define STEVOL	// simple stellar evolution

#define TIMING

#define DTMAXPOWER      -3.0
#define DTMINPOWER     -23.0

/*
 -3.0 	0.125
 -4.0   0.0625
 -5.0   0.03125
 -7.0 	~1e-2
-10.0	~1e-3
.............
-20.0 	~1e-6
-23.0	~1e-7
-26.0	~1e-8
-30.0	~1e-9
*/

//#define DEF_DC

//#define CMCORR

//#define DEBUG_extra

//#define GRAPESHARINGOLD

/*
#define GRAPESHARING
#define DT_GRAPE_LIMIT 300		// DT_GRAPE in seconds
*/

//#define AMD

/****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <sys/time.h>
#include <sys/resource.h>

//#include <grape6.h>
#include "grape6.h"

/* Some "good" functions and constants... */
#define SIG(x)   ( ((x)<0) ? (-1):(1) )
#define ABS(x)   ( ((x)<0) ? (-x):(x) )
#define MAX(a,b) ( ((a)>(b)) ? (a):(b) )
#define MIN(a,b) ( ((a)<(b)) ? (a):(b) )
#define SQR(x)   ( (x)*(x) )

#define Pi      3.141592653589793238462643

#ifdef NORM
#define G     	6.6704E-11			// SI
#define R_gas 	8.31434E+03			// SI
#define k_gas 	1.38066E-23			// SI
#define Mp    	1.6726E-27			// kg
#define Msol  	1.9891E+30			// kg
#define AU  	1.49597870691E+11		// m
#define pc    	3.085678E+16			// m
#define kpc   	(1.0E+03*3.085678E+16)		// m
#define km    	1.0E+03				// m
#define cm3   	1.0E-06				// m^3
#define Year  	3.1556926E+07			// sec
#define Myr  	(1.0E+06*3.1556926E+07)		// sec
#endif

/*
    1KB =	1024
    2KB = 	2048
    4KB = 	4096
    8KB = 	8192
   16KB =      16384
   32KB =      32768
   64KB =      65536
  128KB =     131072
  256KB =     262144
  512KB =     524288
 1024KB =    1048576   ->   1MB
*/

#define KB         1024
#define N_MAX      128*KB

int    i, j, jjj, k, ni, nj, iii, 
       diskstep=0, power, 
       tmp_i;

double dt_disk, dt_contr,
       t_disk=0.0, t_contr=0.0, min_t=0.0,
       t_end, time_cur, dt_min, dt_max, dt_new,
       eta_s, eta,
       E_pot, E_kin, E_pot_ext, 
       E_tot, E_tot_0, DE_tot, 
       A_drag,
       mcm, rcm_mod, vcm_mod,
       rcm_sum=0.0, vcm_sum=0.0,
       eps=0.0, eps2,
       h2, xcm[3], vcm[3], mom[3],
       a2_mod, adot2_mod,
       dt_tmp, dt2half, dt3over6, dt4over24, dt5over120,
       dtinv, dt2inv, dt3inv,
       a0mia1, ad04plad12, ad0plad1,
       a2[3], a3[3], a2dot1[3],
       a1abs, adot1abs, a2dot1abs, a3dot1abs,
       Timesteps=0.0, n_act_sum=0.0, n_act_distr[N_MAX], g6_calls=0.0,
       over2=1.0/2.0, over3=1.0/3.0, over6=1.0/6.0,
       tmp1, tmp2, tmp;

char   inp_fname[30], out_fname[30];


/* global variables */

int    N, ind[N_MAX];
double m[N_MAX], x[N_MAX][3], v[N_MAX][3],
       pot[N_MAX], a[N_MAX][3], adot[N_MAX][3],
       t[N_MAX], dt[N_MAX];


/* data for active particles */

int    n_act, ind_act[N_MAX];
double m_act[N_MAX],
       x_act[N_MAX][3], v_act[N_MAX][3],
       pot_act[N_MAX], a_act[N_MAX][3], adot_act[N_MAX][3],
       t_act[N_MAX], dt_act[N_MAX],
       x_act_new[N_MAX][3], v_act_new[N_MAX][3],
       pot_act_new[N_MAX], a_act_new[N_MAX][3], adot_act_new[N_MAX][3];

double pot_act_tmp[N_MAX], a_act_tmp[N_MAX][3], adot_act_tmp[N_MAX][3];

FILE   *inp, *out, *tmp_file, *dbg;


/* for simple timing... */

double CPU_time_real0, CPU_time_user0, CPU_time_syst0;
double CPU_time_real,  CPU_time_user,  CPU_time_syst;


/* for extended timing... */

#ifdef TIMING

double CPU_tmp_real0, CPU_tmp_user0, CPU_tmp_syst0;
double CPU_tmp_real,  CPU_tmp_user,  CPU_tmp_syst;

double DT_TOT, 
       DT_ACT_DEF, DT_ACT_PRED, 
       DT_ACT_GRAV, DT_EXT_GRAV,
       DT_ACT_CORR, DT_ACT_LOAD, 
       DT_STEVOL, DT_STARDISK, DT_STARDESTR;

#endif


#ifdef DEF_DC

int    l;
double R_LIM2; 
double R_LIM[] = { 1000, 500, 100, 50, 10, 5, 1, 0.5 };
double x_tmp[N_MAX][3], v_tmp[N_MAX][3];

#endif

double mcm_tmp, xcm_tmp[3], vcm_tmp[3];


double GRAPE_time0, GRAPE_time, dt_grape;


/* some settings for GRAPE board */ 

int clusterid=0, npipe=48, 
    new_tunit=51, new_xunit=51,
    ii, nn, index_i[48];

int aflag=1, jflag=1, pflag=1;

double ti=0.0,
       a2by18[3], a1by6[3], aby2[3],
       h2_i[48], x_i[48][3], v_i[48][3],
       p_i[48], a_i[48][3], jerk_i[48][3];


/* normalization... */
       
#ifdef NORM
double m_norm, r_norm, v_norm, t_norm;
#endif


/* external potential... */

//#ifdef EXTPOT
/*
double m_bulge, a_bulge, b_bulge,
       m_disk,	a_disk,	b_disk,
       m_halo,	a_halo,	b_halo,
       x_ext,	y_ext,	z_ext,
       vx_ext,	vy_ext,	vz_ext, 
       x_ij,  y_ij,  z_ij,  
       x2_ij, y2_ij, z2_ij, 
       vx_ij, vy_ij, vz_ij, rv_ij,     
       r_tmp, r2_tmp, z_tmp, z2_tmp;
*/

int    num_bh = 0;

double m_bh, b_bh, e_bh = 0.0, 
       r, r2, 
       x_ij,  y_ij,  z_ij,  
       vx_ij, vy_ij, vz_ij, rv_ij;
       
double pot_ext[N_MAX], a_ext[N_MAX][3], adot_ext[N_MAX][3];
       
//#endif


/* stellar evolution... */

//#ifdef STEVOL
int    num_dead, event[N_MAX];
double dt_stevol, t_stevol;
double m0[N_MAX], ZZZ[N_MAX], t0[N_MAX];
//#endif


/* STARDISK Drag force... */

#ifdef STARDISK

double /*pot_drag[N_MAX],*/ a_drag[N_MAX][3], adot_drag[N_MAX][3];

double Q_DRAG, HZ, D_MASS, D_DELTA; // parameters of the gas-disk
double z_new_drag;
double TWOPI;
double D_SCALE;	/* Inner cut-off (lagr radius for 10%) dlja m_bh = 0.1 */
double D_SIGM;

#include "drag_force.c"

#endif


#ifdef DEBUG

int    ind_sort[N_MAX];

double dt_sort[N_MAX];

/***************************************************************/
void d_swap(double *a, double *b)
{
register double tmp;
tmp = *a; *a = *b; *b = tmp;
}

void i_swap(int *a, int *b)
{
register int tmp;
tmp = *a; *a = *b; *b = tmp;
}
/***************************************************************/
/***************************************************************/
void my_sort(int l, int r, double *arr, int *ind)
{

int    i, j, cikl;
double tmp;

i = l; j = r;
tmp = arr[(l+r)/2];

cikl = 1;

while (cikl)
  {
    while (arr[i]<tmp) i++;
    while (tmp<arr[j]) j--;

    if (i<=j)
     {
	d_swap(&arr[i],&arr[j]);
	i_swap(&ind[i],&ind[j]);
	i++; j--;
     }
    else
     {
	cikl = 0;
     }

  }
				      
  if (l<j) my_sort(l, j, arr, ind);
  if (i<r) my_sort(i, r, arr, ind);				      
}
/***************************************************************/
/***************************************************************/
double my_select(int n_1, int n_2, int k, double *arr, int *ind)
{

int   i, ir, j, l, mid, a_ind;
double a;

l  = n_1;
ir = n_2;

for(;;)
  {

  if (ir <= l+1) 
    {

    if (ir == l+1 && arr[ir] < arr[l]) 
      {
      d_swap(&arr[l],&arr[ir]);
      i_swap(&ind[l],&ind[ir]);
      }

    return(arr[k]);

    } 
    else 
    {

    mid=(l+ir) >> 1;

    d_swap(&arr[mid],&arr[l+1]);
    i_swap(&ind[mid],&ind[l+1]);

    if (arr[l+1] > arr[ir]) 
      {
      d_swap(&arr[l+1],&arr[ir]);
      i_swap(&ind[l+1],&ind[ir]);
      }

    if (arr[l] > arr[ir]) 
      {
      d_swap(&arr[l],&arr[ir]);
      i_swap(&ind[l],&ind[ir]);
      }

    if (arr[l+1] > arr[l])
      {
      d_swap(&arr[l+1],&arr[l]);
      i_swap(&ind[l+1],&ind[l]);
      }

    i = l+1;
    j = ir;
    a = arr[l];
    a_ind = ind[l];

    for (;;) 
     {
     do i++; while (arr[i] < a);
     do j--; while (arr[j] > a);
     if (j < i) break;
     d_swap(&arr[i],&arr[j]);
     i_swap(&ind[i],&ind[j]);
     }

    arr[l] = arr[j];
    ind[l] = ind[j];
    arr[j] = a;
    ind[j] = a_ind;
    if (j >= k) ir = j-1;
    if (j <= k) l = i;

    }
  }
}
/***************************************************************/
#endif

/****************************************************************************/
/****************************************************************************/
/* RAND_MAX = 2147483647 */
/* my_rand  :  0.0 - 1.0 */
/* my_rand2 : -1.0 - 1.0 */

double my_rand(void)
{
return( (double)(rand()/(double)RAND_MAX) );
}

double my_rand2(void)
{
return (double)(2.0)*((rand() - RAND_MAX/2)/(double)RAND_MAX);
}
/****************************************************************************/

#ifdef STARDESTR
#include "star_destr.c"
#endif
       
/****************************************************************************/
void get_CPU_time(double *time_real, double *time_user, double *time_syst)
  {
  struct rusage xxx;
  double sec_u, microsec_u, sec_s, microsec_s;

  getrusage(RUSAGE_SELF,&xxx);

  sec_u = xxx.ru_utime.tv_sec;
  sec_s = xxx.ru_stime.tv_sec;

  microsec_u = xxx.ru_utime.tv_usec;
  microsec_s = xxx.ru_stime.tv_usec;

  *time_user = sec_u + microsec_u * 1.0E-06;
  *time_syst = sec_s + microsec_s * 1.0E-06; 

  *time_real = time(NULL);
  }
/****************************************************************************/
/****************************************************************************/
double mod( double x, double y )
  {
  return( ((x/y)-((int)(x/y)))/y );
  }
/****************************************************************************/

#ifdef STEVOL
/****************************************************************************/
double t_dead(double m, double Z)
/*
m - mass of the star (in Msol)
Z - metallicity (e.g. Zsol = 0.02)
t - return value of star lifetime (in Year)

Raiteri, Villata & Navarro, 1996, A&A, 315, 105
*/
{
  double a0, a1, a2, lZ, lm, tmp=0.0;
  
  m *= (m_norm/Msol);

  lZ = log10(Z);
  lm = log10(m);

  a0 = 10.130 + 0.07547*lZ - 0.008084*lZ*lZ;
  a1 = -4.424 - 0.79390*lZ - 0.118700*lZ*lZ;
  a2 =  1.262 + 0.33850*lZ + 0.054170*lZ*lZ;

  tmp = a0 + a1*lm + a2*lm*lm;

  tmp = pow(10.0,tmp);

  tmp *= (Year/t_norm);

  return(tmp);
}
/****************************************************************************/

/****************************************************************************/
double m_loss(double m, double Z)
/*
m     - initial mass of the star (in Msol)
Z     - metallicity (e.g. Zsol = 0.02)
m_ret - returned maass to ISM from the star (in Msol)

approx. from data of van den Hoek & Groenewegen, 1997, A&AS, 123, 305
*/
{
  double tmp=0.0;
  
  m *= (m_norm/Msol);

  tmp = (-0.46168 + 0.912274 * m);	// first approx.

//  tmp = -0.420542*pow(Z, -0.0176601) + (0.901495 + 0.629441*Z) * m;

  tmp *= (Msol/m_norm);

  return(tmp);
}
/****************************************************************************/

/****************************************************************************/
double m_curr(double m0, double Z, double t0, double t)
{
  double td, ml, tmp=0.0;
  
  td = t_dead(m0, Z);
  ml = m_loss(m0, Z);
  
  if( (t-t0) < td )
    {
//    tmp = m0;					// instant mass loss
    tmp = m0 - (t-t0)/(td-t0) * ml;		// linear mass loss
    }
  else
    {
    tmp = m0 - ml;
    
    if(event[i] == 0) 
      {
      num_dead++;
      event[i] = 1;
      }
      
    } /* if( (t-t0) < td ) */

  return(tmp);
}
/****************************************************************************/
#endif

#ifdef DEF_DC
/****************************************************************************/
void def_dc(int N,         double m[],
            double x[][3], double v[][3],
	    double xcm[],  double vcm[])
{

  for(i=0; i<N; i++)
    {
  for(k=0; k<3; k++)
    {
    x_tmp[i][k] = x[i][k];
    v_tmp[i][k] = v[i][k];
    } /* k */
    } /* i */


  for(k=0; k<3; k++)
    {
    xcm[k] = 0.0;
    vcm[k] = 0.0;
    } /* k */  


  for(l=0; l<sizeof(R_LIM)/sizeof(R_LIM[0]); l++)
    {

    R_LIM2 = SQR(R_LIM[l]);

    mcm_tmp = 0.0;

    for(k=0;k<3;k++) 
      {
      xcm_tmp[k] = 0.0; 
      vcm_tmp[k] = 0.0;
      } /* k */

    for(i=0; i<N; i++)
      {
      
      tmp = SQR(x_tmp[i][0]) + SQR(x_tmp[i][1]) + SQR(x_tmp[i][2]);

      if(tmp < R_LIM2)
        {
      	
        mcm_tmp += m[i];
      
        for(k=0;k<3;k++) 
          {
          xcm_tmp[k] += m[i] * x_tmp[i][k]; 
          vcm_tmp[k] += m[i] * v_tmp[i][k];
          } /* k */
        
        } /* if(tmp < R_LIM2) */
        
      }  /* i */


    for(k=0;k<3;k++) 
      {
      xcm_tmp[k] /= mcm_tmp; 
      vcm_tmp[k] /= mcm_tmp;
      } /* k */


    for(i=0; i<N; i++)
      {        	
    for(k=0;k<3;k++) 
      {
      x_tmp[i][k] -= xcm_tmp[k]; 
      v_tmp[i][k] -= vcm_tmp[k];
      } /* k */  
      } /* i */


    for(k=0;k<3;k++) 
      {
      xcm[k] += xcm_tmp[k];
      vcm[k] += vcm_tmp[k];
      } /* k */
    
    
    } /* Do sizeof(R_LIM)/sizeof(R_LIM[0]) iterations for R_LIM */
    
}

/****************************************************************************/
#endif

/****************************************************************************/
void read_data()
{
  inp = fopen(inp_fname,"r");
  fscanf(inp,"%d \n",  &diskstep);
  fscanf(inp,"%d \n",  &N);
  fscanf(inp,"%lE \n", &time_cur);
  for(i=0; i<N; i++) 
  fscanf(inp,"%d %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %d \n", 
             &ind[i], &m[i], &x[i][0], &x[i][1], &x[i][2], &v[i][0], &v[i][1], &v[i][2], &m0[i], &ZZZ[i], &t0[i], &event[i]);
  fclose(inp);  
}
/****************************************************************************/

/****************************************************************************/
void cm_corr()
{    
    mcm = 0.0;

    for(k=0;k<3;k++) 
      {
      xcm[k] = 0.0; 
      vcm[k] = 0.0;
      } /* k */
	
    for(i=0; i<N; i++)
      {
        	
      mcm += m[i];
        
      for(k=0;k<3;k++) 
        {
        xcm[k] += m[i] * x[i][k]; 
        vcm[k] += m[i] * v[i][k];
        } /* k */
          
      }  /* i */

    for(k=0;k<3;k++) 
      {
      xcm[k] /= mcm; 
      vcm[k] /= mcm;
      } /* k */

    for(i=0; i<N; i++)
      {
        	
      for(k=0;k<3;k++) 
        {
        x[i][k] -= xcm[k]; 
        v[i][k] -= vcm[k];
        } /* k */
        
      } /* i */
}
/****************************************************************************/

/****************************************************************************/
void write_snap_data()
{
  sprintf(out_fname,"%04d.dat",diskstep);
  out = fopen(out_fname,"w");
  fprintf(out,"%04d \n", diskstep);
  fprintf(out,"%06d \n", N);
  fprintf(out,"%.6E \n", time_cur);
  for(i=0; i<N; i++) 
  fprintf(out,"%06d  %.6E  % .6E % .6E % .6E  % .6E % .6E % .6E  %.6E %.6E %.6E %d \n", 
  	       ind[i], 
  	       m[i], 
  	       x[i][0], x[i][1], x[i][2], 
  	       v[i][0], v[i][1], v[i][2], 
  	       m0[i], ZZZ[i], t0[i], event[i]);
  fclose(out);    
}
/****************************************************************************/

/****************************************************************************/
void write_cont_data()
{
  out = fopen("data.con","w");
  fprintf(out,"%04d \n", diskstep);
  fprintf(out,"%06d \n", N);
  fprintf(out,"%.9E \n", time_cur);
  for(i=0; i<N; i++) 
  fprintf(out,"%06d  %.9E  % .9E % .9E % .9E  % .9E % .9E % .9E  %.9E %.9E %.9E %d \n", 
  	       ind[i], 
  	       m[i], 
  	       x[i][0], x[i][1], x[i][2], 
  	       v[i][0], v[i][1], v[i][2], 
  	       m0[i], ZZZ[i], t0[i], event[i]);
  fclose(out);
}
/****************************************************************************/


#ifdef DEBUG

/****************************************************************************/
void write_dbg_data()
{

  double tmp_r, tmp_v, tmp_a, tmp_adot, tmp_a_ext, tmp_adot_ext, tmp_a_drag, tmp_adot_drag;
  
  for(i=0; i<N; i++)
    {
    ind_sort[i] = ind[i];
    dt_sort[i] = dt[i];
    }
  
  my_sort(0, N-1, dt_sort, ind_sort);
   
  //tmp = my_select(0, N-1, Nb-1, d, ind); d[Nb-1] = tmp;

  sprintf(out_fname,"%04d.dbg",diskstep);

  out = fopen(out_fname,"w");

  fprintf(out,"%04d \n", diskstep);
  fprintf(out,"%06d \t %06d \n", N, n_act);
  fprintf(out,"%.6E \n", time_cur);

  for(j=0; j<N; j++) 
    {
    i = ind_sort[j];

    tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) );
    tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) );

    tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][1]) );
    tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][1]) );

    tmp_a_ext = sqrt( SQR(a_ext[i][0]) + SQR(a_ext[i][1]) + SQR(a_ext[i][1]) );
    tmp_adot_ext = sqrt( SQR(adot_ext[i][0]) + SQR(adot_ext[i][1]) + SQR(adot_ext[i][1]) );

    tmp_a_drag = sqrt( SQR(a_drag[ind_act[i]][0]) + SQR(a_drag[ind_act[i]][1]) + SQR(a_drag[ind_act[i]][1]) );
    tmp_adot_drag = sqrt( SQR(adot_drag[ind_act[i]][0]) + SQR(adot_drag[ind_act[i]][1]) + SQR(adot_drag[ind_act[i]][1]) );
    
    fprintf(out,"%06d %06d %06d \t %.6E \t % .6E % .6E % .6E %.6E \t % .6E % .6E % .6E %.6E \t % .6E \t % .6E % .6E % .6E %.6E \t % .6E % .6E % .6E %.6E \t % .6E \t % .6E % .6E % .6E %.6E \t % .6E % .6E % .6E %.6E \t % .6E % .6E % .6E %.6E \t % .6E % .6E % .6E %.6E \t %.6E \n", 
  		 ind[i], ind_act[i], ind_sort[j], 
  		 m[i], 
  		 x[i][0], x[i][1], x[i][2], tmp_r,
  		 v[i][0], v[i][1], v[i][2], tmp_v,
		 pot[i], 
  		 a[i][0], a[i][1], a[i][2], tmp_a,
		 adot[i][0], adot[i][1], adot[i][2], tmp_adot, 
		 pot_ext[i], 
  		 a_ext[i][0], a_ext[i][1], a_ext[i][2], tmp_a_ext,
		 adot_ext[i][0], adot_ext[i][1], adot_ext[i][2], tmp_adot_ext, 
		 //pot_drag[i], 
  		 a_drag[ind_act[i]][0], a_drag[ind_act[i]][1], a_drag[ind_act[i]][2], tmp_a_drag,
		 adot_drag[ind_act[i]][0], adot_drag[ind_act[i]][1], adot_drag[ind_act[i]][2], tmp_adot_drag, 
		 dt_sort[j]);
		
    }

  fclose(out);    

}
/****************************************************************************/

#endif

/****************************************************************************/
void calc_self_grav()
{

  for(i=0; i<ni; i+=npipe)
    {

    nn = npipe;
    if(ni-i < npipe) nn = ni - i;

    for(ii=0; ii<nn; ii++)
      {
      iii = ind_act[i+ii];

      index_i[ii] = ind_act[i+ii];
      h2_i[ii]    = h2;

      for(k=0;k<3;k++)
        {
        x_i[ii][k] = x_act_new[i+ii][k]-xcm[k];
        v_i[ii][k] = v_act_new[i+ii][k]-vcm[k];
        } /* k */

      p_i[ii] = pot_act_tmp[iii];

      for(k=0;k<3;k++)
        {
        a_i[ii][k] = a_act_tmp[iii][k];
        jerk_i[ii][k] = adot_act_tmp[iii][k];
        } /* k */

      } /* ii */

      g6calc_firsthalf(clusterid, N, nn, index_i, x_i, v_i , a_i, jerk_i, p_i, eps2, h2_i);
      g6calc_lasthalf(clusterid, N, nn, index_i, x_i, v_i, eps2, h2_i, a_i, jerk_i, p_i);

      g6_calls++;


      if(Timesteps == 0.0)
        {
        for(ii=0; ii<nn; ii++)
          {
          if(ABS(jerk_i[ii][0]) < 1.0E-05) jerk_i[ii][0] = 1.0E-05;
          if(ABS(jerk_i[ii][1]) < 1.0E-05) jerk_i[ii][1] = 1.0E-05;
          if(ABS(jerk_i[ii][2]) < 1.0E-05) jerk_i[ii][2] = 1.0E-05;
	
          if(ABS(jerk_i[ii][0]) > 100.0) jerk_i[ii][0] = 100.0;
          if(ABS(jerk_i[ii][1]) > 100.0) jerk_i[ii][1] = 100.0;
          if(ABS(jerk_i[ii][2]) > 100.0) jerk_i[ii][2] = 100.0;
	  }
	
        g6calc_firsthalf(clusterid, N, nn, index_i, x_i, v_i , a_i, jerk_i, p_i, eps2, h2_i);
        g6calc_lasthalf(clusterid, N, nn, index_i, x_i, v_i, eps2, h2_i, a_i, jerk_i, p_i);
      
        g6_calls++;
        }


      for(ii=0; ii<nn; ii++) 
        {
        pot_act_new[i+ii] = p_i[ii];

        for(k=0;k<3;k++) 
          {
          a_act_new[i+ii][k] = a_i[ii][k];
          adot_act_new[i+ii][k] = jerk_i[ii][k];
          } /* k */
          
        } /* ii */

      } /* i */


  /* Store the new value of the local force etc... */

  for(i=0; i<n_act; i++)
    {
    iii = ind_act[i];
    
    pot_act_tmp[iii] = pot_act_new[i];

    for(k=0;k<3;k++) 
      {
      a_act_tmp[iii][k] = a_act_new[i][k];
      adot_act_tmp[iii][k] = adot_act_new[i][k];
      } /* k */

    } /* i */

}
/****************************************************************************/

#ifdef EXTPOT
/****************************************************************************/
void calc_ext_grav()
{
	
  for(i=0; i<n_act; i++) {
    
    pot_ext[i] = 0.0;

    for(k=0;k<3;k++)   {
      a_ext[i][k] = 0.0;
      adot_ext[i][k] = 0.0;
    }	
    
  }


  for(i=0; i<n_act; i++)
    {

    /* For simple Plummer potential... */

    x_ij = x_act_new[i][0];
    y_ij = x_act_new[i][1];
    z_ij = x_act_new[i][2];
  
    vx_ij = v_act_new[i][0];
    vy_ij = v_act_new[i][1];
    vz_ij = v_act_new[i][2];

    r2 = SQR(b_bh);
        
    r2 += x_ij*x_ij + y_ij*y_ij + z_ij*z_ij;   r = sqrt(r2);
      
    rv_ij = vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij;

    tmp = m_bh / r;

    pot_ext[i] -= tmp;

    tmp /= r2;

    a_ext[i][0] -= tmp * x_ij;
    a_ext[i][1] -= tmp * y_ij;
    a_ext[i][2] -= tmp * z_ij;

    adot_ext[i][0] -= tmp * (vx_ij - 3.0*rv_ij * x_ij/r2);
    adot_ext[i][1] -= tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
    adot_ext[i][2] -= tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);

    } /* i */
    
}
/****************************************************************************/
#endif		// EXTPOT

/****************************************************************************/
void energy_contr()
{

  E_pot = 0.0;
  for(i=0; i<N; i++) E_pot += m[i]*pot[i];
  E_pot *= 0.5;
  
  E_pot_ext = 0.0;
  for(i=0; i<N; i++) E_pot_ext += m[i]*pot_ext[i];

  E_kin = 0.0;
  for(i=0; i<N; i++) E_kin += m[i]*( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
  E_kin *= 0.5;


  if(Timesteps == 0.0)
    {
    rcm_sum = 0.0; 
    vcm_sum = 0.0;
    }

  
  mcm     = 0.0; 
  rcm_mod = 0.0; 
  vcm_mod = 0.0; 

  for(k=0;k<3;k++) 
    {
    xcm[k] = 0.0; 
    vcm[k] = 0.0;
    } /* k */

  for(i=0; i<N; i++)
    {
    mcm += m[i];

    for(k=0;k<3;k++) 
      {
      xcm[k] += m[i] * x[i][k]; 
      vcm[k] += m[i] * v[i][k];
      } /* k */
        
    } /* i */

  for(k=0;k<3;k++) 
    {
    xcm[k] /= mcm; 
    vcm[k] /= mcm;
    } /* k */


  rcm_mod = sqrt(xcm[0]*xcm[0]+xcm[1]*xcm[1]+xcm[2]*xcm[2]);
  vcm_mod = sqrt(vcm[0]*vcm[0]+vcm[1]*vcm[1]+vcm[2]*vcm[2]);
  
  rcm_sum += rcm_mod;
  vcm_sum += vcm_mod;
  

  for(k=0;k<3;k++) mom[k] = 0.0;

  for(i=0; i<N; i++)
    {      	
    mom[0] += m[i] * ( x[i][1] * v[i][2] - x[i][2] * v[i][1] );
    mom[1] += m[i] * ( x[i][2] * v[i][0] - x[i][0] * v[i][2] );
    mom[2] += m[i] * ( x[i][0] * v[i][1] - x[i][1] * v[i][0] );
    } /* i */


  get_CPU_time(&CPU_time_real, &CPU_time_user, &CPU_time_syst);

  E_tot = E_pot + E_pot_ext + E_kin + e_bh + A_drag;

  if(Timesteps == 0.0) E_tot_0 = E_tot;

  DE_tot = E_tot - E_tot_0;

  printf("%.4E   %.3E %.3E %.3E   % .4E % .4E % .4E % .4E % .4E % .6E   % .4E   %.4E \n",
         time_cur, Timesteps, n_act_sum, g6_calls,
         E_pot, E_pot_ext, E_kin, e_bh, A_drag, E_tot, DE_tot, 
         CPU_time_user-CPU_time_user0);


  fflush(stdout);

  out = fopen("contr.dat","a");


  fprintf(out,"%.8E  %.8E  %.8E  %.8E   % .8E % .8E % .8E % .8E % .8E % .8E % .8E   % .8E % .8E   % .8E % .8E % .8E   %.8E %.8E %.8E \n",
          time_cur, Timesteps, n_act_sum, g6_calls,
          E_pot, E_pot_ext, E_kin, e_bh, A_drag, E_tot, DE_tot,
//          rcm_mod, vcm_mod,
          rcm_sum, vcm_sum,
          mom[0], mom[1], mom[2],
          CPU_time_real-CPU_time_real0, CPU_time_user-CPU_time_user0, CPU_time_syst-CPU_time_syst0);

    
  fclose(out);
  
  E_tot_0 = E_tot;

}
/****************************************************************************/

#ifdef GRAPESHARING
/****************************************************************************/
void grape_sharing()
{
    GRAPE_time = time(NULL);
  
    dt_grape = GRAPE_time - GRAPE_time0;

    if(dt_grape >= DT_GRAPE_LIMIT)
      {

      /* close the GRAPE */

      g6_close(clusterid);

      /* Give the chance to catch the GRAPE to someone else... :-) */

      /* init the GRAPE */

      g6_open(clusterid);   GRAPE_time0 = time(NULL);

      npipe = g6_npipes();

      g6_set_tunit(new_tunit);
      g6_set_xunit(new_xunit);

      g6_set_overflow_flag_test_mode(aflag, jflag, pflag);

#ifdef AMD
      g6_initialize_jp_buffer(clusterid, 10000);
#endif

      /* load the new values for particles to the GRAPE */

      nj=N; 

      /* load the nj particles to the GRAPE */

      for(j=0; j<nj; j++) 
        {
        for(k=0;k<3;k++) 
          {
          a2by18[k] = 0.0; a1by6[k] = adot[j][k]*over6; aby2[k] = a[j][k]*over2;
          } /* k */

        g6_set_j_particle(clusterid, j, ind[j], t[j], dt[j], m[j], a2by18, a1by6, aby2, v[j], x[j]);
        } /* j */

#ifdef AMD
    g6_flush_jp_buffer(clusterid);
#endif

      } /* if(dt_grape >= DT_GRAPE_LIMIT) */
}
/****************************************************************************/
#endif

#ifdef GRAPESHARINGOLD
/****************************************************************************/
void grape_sharing_old()
{
    /* close the GRAPE */

    g6_close(clusterid);

    /* Give the chance to catch the GRAPE to someone else... :-) */

    /* init the GRAPE */
   
    g6_open(clusterid);

    npipe = g6_npipes();

    g6_set_tunit(new_tunit);
    g6_set_xunit(new_xunit);

    g6_set_overflow_flag_test_mode(aflag, jflag, pflag);

#ifdef AMD
    g6_initialize_jp_buffer(clusterid, 10000);
#endif

    /* load the new values for particles to the GRAPE */

    nj=N; 

    /* load the nj particles to the GRAPE */

    for(j=0; j<nj; j++) 
      {
      for(k=0;k<3;k++) 
        {
        a2by18[k] = 0.0; a1by6[k] = adot[j][k]*over6; aby2[k] = a[j][k]*over2;
        } /* k */

      g6_set_j_particle(clusterid, j, ind[j], t[j], dt[j], m[j], a2by18, a1by6, aby2, v[j], x[j]);
      } /* j */

#ifdef AMD
    g6_flush_jp_buffer(clusterid);
#endif
}
/****************************************************************************/
#endif


/****************************************************************************/
/****************************************************************************/
/****************************************************************************/

int main(int argc, char *argv[])
  {

  /* INIT the rand() !!! */
  srand(19640916);                 /* it is just my birthday :-) */
  /* srand(time(NULL)); */ 


  /* Read the basic input parameters of the run */

  inp = fopen("phi-GRAPE.cfg","r");
  fscanf(inp,"%lE %lE %lE %lE %lE %s", &eps, &t_end, &dt_disk, &dt_contr, &eta, inp_fname);
  fclose(inp);

  /*
  eps		: Plummer softening parameter (can be even 0)
  t_end      	: end time of calculation
  dt_disk  	: interval of snapshot files output (0xxx.dat)
  dt_contr 	: interval for the energy control output (contr.dat)
  eta	     	: parameter for timestep determination
  inp_data 	: name of the input file (data.inp)
  */


   /* Input of the gas-disc parameters  */
#ifdef STARDISK
  inp = fopen("drag.ini","r");
  fscanf(inp,"%lE %lE %lE %lE", &Q_DRAG, &HZ, &D_MASS, &D_DELTA);
  fclose(inp);

  TWOPI = 2.0*Pi;
  D_SCALE = 0.31;	/* Inner cut-off (lagr radius for 10%) dlja m_bh = 0.1 */
  D_SIGM = 5.0*D_MASS/(4*TWOPI*sqrt(TWOPI)*HZ*SQR(D_SCALE));
#endif


  /* Normalization... */
#ifdef NORM
  inp = fopen("phi-GRAPE.norm","r");
  fscanf(inp,"%lE %lE", &m_norm, &r_norm);
  fclose(inp);
  
  m_norm *= Msol;
  r_norm *= pc;

  v_norm = sqrt(G*m_norm/r_norm);
  t_norm = (r_norm/v_norm);
#endif		// NORM


  /* Read the data for EXTernal POTential */
  
#ifdef EXTPOT
  inp = fopen("phi-GRAPE.ext","r");
  fscanf(inp,"%lE %lE", &m_bh, &b_bh);
  fclose(inp);
#endif		// EXTPOT



  /* read the global data for particles */
  
  read_data();


  /* possible CM correction of the initial datafile... */
#ifdef CMCORR

  if(diskstep == 0) cm_corr();

#endif		// CMCORR

        
  printf("\n");
  printf("Begin the calculation of phi-GRAPE (BH+STARDISK) program (serial) \n"); 
  printf("\n");
  printf("N       = %06d \t eps      = %.6E \n", N, eps);
  printf("t_beg   = %.6E \t t_end    = %.6E \n", time_cur, t_end);
  printf("dt_disk = %.6E \t dt_contr = %.6E \n", dt_disk, dt_contr);
  printf("eta     = %.6E \n", eta);
  printf("\n"); 

#ifdef NORM
  printf("Normalization: \n"); 
  printf("\n");
  printf("m_norm = %.4E [Msol]  r_norm = %.4E [pc] \n", m_norm/Msol, r_norm/pc);
  printf("v_morm = %.4E [km/s]  t_morm = %.4E [Myr] \n", v_norm/km, t_norm/Myr);
  printf("\n");  
#endif		// NORM

#ifdef EXTPOT
  printf("External Potential: \n"); 
  printf("\n");
  printf("m_bh    = %.6E \t b_bh     = %.6E \n", m_bh, b_bh);
  printf("\n");
#endif		// EXTPOT

  fflush(stdout);



  if(diskstep == 0)
    {

    write_snap_data();

    write_cont_data();

    out = fopen("contr.dat","w");
    fclose(out);

#ifdef TIMING
    out = fopen("timing.dat","w");
    fclose(out);
#endif		// TIMING

#ifdef STARDESTR
    out = fopen("mass-bh.dat","w");
    fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh);
    fclose(out);
#endif		// STARDESTR
    
    }




  get_CPU_time(&CPU_time_real0, &CPU_time_user0, &CPU_time_syst0);


//  eta_s = eta/32.0; 
  eta_s = eta/32.0;

  eps2 = SQR(eps);
  h2 = eps2;

  dt_min = 1.0*pow(2.0, DTMINPOWER);
  dt_max = 1.0*pow(2.0, DTMAXPOWER);

  if(dt_disk == dt_contr)
    dt_max = dt_contr;
  else
    dt_max = MIN(dt_disk, dt_contr); 


#ifdef STEVOL    
  dt_stevol = dt_max;
  t_stevol  = time_cur + dt_stevol;  
#endif		// STEVOL


  t_disk  = time_cur + dt_disk;
  t_contr = time_cur + dt_contr;


  for(i=0; i<N; i++) 
    {
    t[i]  = time_cur;
    dt[i] = dt_min;
    }



  /* init the GRAPE */

  g6_open(clusterid);   GRAPE_time0 = time(NULL);

  npipe = g6_npipes();

  g6_set_tunit(new_tunit);
  g6_set_xunit(new_xunit);

  g6_set_overflow_flag_test_mode(aflag, jflag, pflag);

#ifdef AMD
  g6_initialize_jp_buffer(clusterid, 10000);
#endif



#ifdef DEF_DC
  def_dc(N, m, x, v, xcm, vcm);
#else
  for(k=0;k<3;k++)
    {
    xcm[k] = 0.0;
    vcm[k] = 0.0;
    } /* k */
#endif		// DEF_DC



  /* initial load the particles to the GRAPE */

  nj=N;

  /* load the nj particles to the GRAPE */

  for(j=0; j<nj; j++) 
    {
    for(k=0;k<3;k++) 
      {
      a2by18[k] = 0.0; a1by6[k] = 0.0; aby2[k] = 0.0; 

      xcm_tmp[k] = x[j][k] - xcm[k]; 
      vcm_tmp[k] = v[j][k] - vcm[k];
      } /* k */

//    g6_set_j_particle(clusterid, j, ind[j], t[j], dt[j], m[j], a2by18, a1by6, aby2, v[j], x[j]);
    g6_set_j_particle(clusterid, j, ind[j], t[j], dt[j], m[j], a2by18, a1by6, aby2, vcm_tmp, xcm_tmp);
    } /* j */

#ifdef AMD
    g6_flush_jp_buffer(clusterid);
#endif




  /* define the all particles as a active for the first time grav calc... */

  n_act = N;

  for(i=0; i<n_act; i++)
    {
      	
    ind_act[i] = ind[i];
    m_act[i]   = m[ind_act[i]];

    for(k=0;k<3;k++) 
      {
      x_act_new[i][k] = x[ind_act[i]][k]; 
      v_act_new[i][k] = v[ind_act[i]][k];
      x_act[i][k] = x[ind_act[i]][k]; 
      v_act[i][k] = v[ind_act[i]][k];
      } /* k */

    t_act[i]  = t[ind_act[i]];
    dt_act[i] = dt[ind_act[i]];
          
    } /* i */


    
  /* Set initial dummy values for pot, a & adot... */

  for(i=0; i<n_act; i++)
    {

    pot_act_tmp[i] = -1.0;

    for(k=0;k<3;k++) 
      {
      a_act_tmp[i][k] = -10.0;
      adot_act_tmp[i][k] = -100.0;
      } /* k */

    } /* i */

    

  /* calc the grav for the active particles */

  g6_set_ti(clusterid, time_cur);
 
  ni = n_act; 



  /* define the phi, a, adot for all the active particles */
  
  calc_self_grav();



#ifdef EXTPOT

  calc_ext_grav();
  
  /* summed up all components... */

  for(i=0; i<n_act; i++)
    {
//    pot_act_new[i] += pot_ext[i];

    for(k=0;k<3;k++) 
      {
      a_act_new[i][k] += a_ext[i][k];
      adot_act_new[i][k] += adot_ext[i][k];
      }	
    }

#endif		// EXTPOT




#ifdef STARDISK

  for(i=0; i<N; i++) for(k=0; k<3; k++) { a_drag[i][k] = 0.0; adot_drag[i][k] = 0.0; }

  A_drag = 0.0;

  calc_drag_force();

//  A_drag = 0.0;

//  for(i=0; i<n_act; i++) A_drag += m_act[i]*pot_drag[i];
  
  /* summed up all components... */

  for(i=0; i<n_act; i++)
    {
    for(k=0;k<3;k++) 
      {
      a_act_new[i][k] += a_drag[ind_act[i]][k];
      adot_act_new[i][k] += adot_drag[ind_act[i]][k];
      }	
    }
    
#endif		// STARDISK




  /* Return back to the global array's... */

  for(i=0; i<n_act; i++)
    {

    pot[ind_act[i]] = pot_act_new[i];

    for(k=0;k<3;k++) 
      {
      a[ind_act[i]][k] = a_act_new[i][k];
      adot[ind_act[i]][k] = adot_act_new[i][k];
      } /* k */
        
    } /* i */




  /* Energy control... */
  
  energy_contr();

 

  /* Define initial timestep for all particles */

  for(i=0; i<N; i++)
    {

    a2_mod = SQR(a[i][0])+SQR(a[i][1])+SQR(a[i][2]);
    adot2_mod = SQR(adot[i][0])+SQR(adot[i][1])+SQR(adot[i][2]);

/*
    a2_mod = a[i][0]*a[i][0]+a[i][1]*a[i][1]+a[i][2]*a[i][2];
    adot2_mod = adot[i][0]*adot[i][0]+adot[i][1]*adot[i][1]+adot[i][2]*adot[i][2];
*/

    if(adot2_mod == 0)
      dt_tmp = eta_s;
    else
      dt_tmp = eta_s*sqrt(a2_mod/adot2_mod);


    power = log(dt_tmp)/log(2.0) - 1;
//    power = log(dt_tmp)/log(2.0);

    dt_tmp = pow(2.0, (double)power);

    if(dt_tmp < dt_min) dt_tmp = dt_min;
    if(dt_tmp > dt_max) dt_tmp = dt_max;

    dt[i] = dt_tmp;      
    } /* i */




#ifdef DEF_DC
  def_dc(N, m, x, v, xcm, vcm);
#else
  for(k=0;k<3;k++)
    {
    xcm[k] = 0.0;
    vcm[k] = 0.0;
    } /* k */
#endif		// DEF_DC


  /* load the new values for particles to the GRAPE */

  nj=N; 

  /* load the nj particles to the GRAPE */

  for(j=0; j<nj; j++) 
    {
    for(k=0;k<3;k++) 
      {
      a2by18[k] = 0.0; a1by6[k] = adot[j][k]*over6; aby2[k] = a[j][k]*over2;

      xcm_tmp[k] = x[j][k] - xcm[k]; 
      vcm_tmp[k] = v[j][k] - vcm[k];
      } /* k */

//    g6_set_j_particle(clusterid, j, ind[j], t[j], dt[j], m[j], a2by18, a1by6, aby2, v[j], x[j]);
    g6_set_j_particle(clusterid, j, ind[j], t[j], dt[j], m[j], a2by18, a1by6, aby2, vcm_tmp, xcm_tmp);
    } /* j */

#ifdef AMD
    g6_flush_jp_buffer(clusterid);
#endif



#ifdef DEBUG
    write_dbg_data();
#endif


  get_CPU_time(&CPU_time_real0, &CPU_time_user0, &CPU_time_syst0); 


  Timesteps = 0.0;
  n_act_sum = 0.0;

  for(i=1; i<N+1; i++) 
    {
    n_act_distr[i-1] = 0.0;
    }

  g6_calls = 0.0;


#ifdef TIMING
  DT_TOT      = 0.0;

  DT_ACT_DEF  = 0.0;
  DT_ACT_PRED = 0.0;
  DT_ACT_GRAV = 0.0;
  DT_EXT_GRAV = 0.0;
  DT_ACT_CORR = 0.0;
  DT_ACT_LOAD = 0.0;

  DT_STEVOL    = 0.0;
  DT_STARDISK  = 0.0;
  DT_STARDESTR = 0.0;
#endif		// TIMING  



#ifdef DEBUG111
  printf("%6.0f \t %06d \t %6.0f \n", Timesteps, n_act, n_act_sum);
  fflush(stdout);
#endif



  /* The main integration loop */

  while(time_cur <= t_end)
//  while( Timesteps < 1e3 )
    {


  /* Define the minimal time and the active particles list */

#ifdef TIMING
  get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

  min_t = t[0]+dt[0];

  for(i=1; i<N; i++)
    {
    if( (t[i]+dt[i]) < min_t ) min_t = t[i]+dt[i];
    } /* i */

  n_act = 0;
  for(i=0; i<N; i++)
    {      	
    if( (t[i]+dt[i]) == min_t ) 
      {
      ind_act[n_act] = ind[i];
      n_act++;
      }        
    } /* i */

  for(i=0; i<n_act; i++)
    {
    m_act[i] = m[ind_act[i]];
      
    for(k=0;k<3;k++) 
      {
      x_act[i][k] = x[ind_act[i]][k];
      v_act[i][k] = v[ind_act[i]][k];
      } /* k */

    t_act[i]  = t[ind_act[i]];
    dt_act[i] = dt[ind_act[i]];

    pot_act[i] = pot[ind_act[i]];

    for(k=0;k<3;k++) 
      {
      a_act[i][k] = a[ind_act[i]][k];
      adot_act[i][k] = adot[ind_act[i]][k];
      } /* k */

    } /* i */

#ifdef TIMING
  get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
  DT_ACT_DEF += (CPU_tmp_user - CPU_tmp_user0);
#endif





  /* make the prediction for the active particles... */

#ifdef TIMING
  get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

  for(i=0; i<n_act; i++)
    {
    dt_tmp = min_t - t_act[i];
    dt2half = over2*dt_tmp*dt_tmp;
    dt3over6 = over3*dt_tmp*dt2half;

    for(k=0;k<3;k++) 
      {
      x_act_new[i][k] = x_act[i][k] + v_act[i][k]*dt_tmp + a_act[i][k]*dt2half + adot_act[i][k]*dt3over6;
      v_act_new[i][k] = v_act[i][k] + a_act[i][k]*dt_tmp + adot_act[i][k]*dt2half;
      } /* k */
    } /* i */

#ifdef TIMING
  get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
  DT_ACT_PRED += (CPU_tmp_user - CPU_tmp_user0);
#endif




  /* calc the new grav for the active particles */

#ifdef TIMING
  get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

  g6_set_ti(clusterid, min_t);
 
  ni = n_act;

  /* define the phi, a, adot for all the active particles */
  
  calc_self_grav();

#ifdef TIMING
  get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
  DT_ACT_GRAV += (CPU_tmp_user - CPU_tmp_user0);
#endif




#ifdef EXTPOT

#ifdef TIMING
  get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

  calc_ext_grav();
  
  /* summed up all components... */

  for(i=0; i<n_act; i++)
    {
//    pot_act_new[i] += pot_ext[i];

    for(k=0;k<3;k++) 
      {
      a_act_new[i][k] += a_ext[i][k];
      adot_act_new[i][k] += adot_ext[i][k];
      }	
    }
  
#ifdef TIMING
  get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
  DT_EXT_GRAV += (CPU_tmp_user - CPU_tmp_user0);
#endif

#endif		// EXTPOT





#ifdef STARDISK

#ifdef TIMING
  get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

  calc_drag_force();  

//  for(i=0; i<n_act; i++) A_drag += m_act[i]*pot_drag[i];
  
  /* summed up all components... */

  for(i=0; i<n_act; i++)
    {
    for(k=0;k<3;k++) 
      {
      a_act_new[i][k] += a_drag[ind_act[i]][k];
      adot_act_new[i][k] += adot_drag[ind_act[i]][k];
      }	
    }
  
#ifdef TIMING
  get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
  DT_STARDISK += (CPU_tmp_user - CPU_tmp_user0);
#endif

#endif		// STARDISK  




  /* correct the active particles positions etc and define the new timesteps for active particles... */

#ifdef TIMING
  get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

  for(i=0; i<n_act; i++)
    {
  
    dt_tmp = min_t - t_act[i];
      
    dt3over6 = dt_tmp*dt_tmp*dt_tmp/6.0;
    dt4over24 = dt3over6*dt_tmp/4.0;
    dt5over120 = dt4over24*dt_tmp/5.0;
      
    dtinv = 1.0/dt_tmp;
    dt2inv = dtinv*dtinv;
    dt3inv = dt2inv*dtinv;
      
    for(k=0; k<3; k++)
      {
      a0mia1 = a_act[i][k]-a_act_new[i][k];
      ad04plad12 = 4.0*adot_act[i][k] + 2.0*adot_act_new[i][k];
      ad0plad1 = adot_act[i][k] + adot_act_new[i][k];
        
      a2[k] = -6.0*a0mia1*dt2inv - ad04plad12*dtinv;
      a3[k] = 12.0*a0mia1*dt3inv + 6.0*ad0plad1*dt2inv;
  
      x_act_new[i][k] += dt4over24*a2[k] + dt5over120*a3[k];
      v_act_new[i][k] += dt3over6*a2[k] + dt4over24*a3[k];
      } /* k */


    a1abs = sqrt(SQR(a_act_new[i][0])+SQR(a_act_new[i][1])+SQR(a_act_new[i][2]));
    adot1abs = sqrt(SQR(adot_act_new[i][0])+SQR(adot_act_new[i][1])+SQR(adot_act_new[i][2]));
      
    for(k=0; k<3; k++) a2dot1[k] = a2[k] + dt_tmp*a3[k];
      
    a2dot1abs = sqrt(SQR(a2dot1[0])+SQR(a2dot1[1])+SQR(a2dot1[2]));
    a3dot1abs = sqrt(SQR(a3[0])+SQR(a3[1])+SQR(a3[2]));
  
    dt_new = sqrt(eta*(a1abs*a2dot1abs+SQR(adot1abs))/(adot1abs*a3dot1abs+SQR(a2dot1abs)));
    
    if(m_act[i] == 0.0) dt_new = dt_contr;

/*  
    a1abs = sqrt(a_act_new[i][0]*a_act_new[i][0]+a_act_new[i][1]*a_act_new[i][1]+a_act_new[i][2]*a_act_new[i][2]);
    adot1abs = sqrt(adot_act_new[i][0]*adot_act_new[i][0]+adot_act_new[i][1]*adot_act_new[i][1]+adot_act_new[i][2]*adot_act_new[i][2]);
      
    for(k=0; k<3; k++) a2dot1[k] = a2[k] + dt_tmp*a3[k];
      
    a2dot1abs = sqrt(a2dot1[0]*a2dot1[0]+a2dot1[1]*a2dot1[1]+a2dot1[2]*a2dot1[2]);
    a3dot1abs = sqrt(a3[0]*a3[0]+a3[1]*a3[1]+a3[2]*a3[2]);

    dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
*/

    if(dt_new < dt_min) dt_tmp = dt_min; 

    if( (dt_new < dt_tmp) && (dt_new > dt_min) )
      {
      power = log(dt_new)/log(2.0) - 1;
//      power = log(dt_new)/log(2.0);

      dt_tmp = pow(2.0, (double)power);
      }	  


#ifdef STARDISK

   // accretion disk criterion
    z_new_drag = x_act_new[i][2] + dt_tmp * v_act_new[i][2];

    if (fabs(x_act_new[i][2]) > HZ*2.2)
      {
      if (z_new_drag * x_act_new[i][2] < 0.0)
	{
        dt_tmp *= 0.125;
        }
      else
        {
        if (fabs(z_new_drag) < HZ*1.8)
	  {
          dt_tmp *= 0.5;
          }
        }
      }

#endif


    if( (dt_new > 2.0*dt_tmp) && 
        (fmod(min_t, 2.0*dt_tmp) == 0.0) && 
        (2.0*dt_tmp <= dt_max) )
      {
      dt_tmp *= 2.0;
      }


    dt_act[i] = dt_tmp;
  
    t_act[i] = min_t;

    pot_act[i] = pot_act_new[i];
      
    for(k=0; k<3; k++)
      { 
         x_act[i][k] =    x_act_new[i][k];
         v_act[i][k] =    v_act_new[i][k];
         a_act[i][k] =    a_act_new[i][k];
      adot_act[i][k] = adot_act_new[i][k];
      } /* k */


    } /* i */



  /* Return back the new coordinates + etc... of active part to the global array's... */

  for(i=0; i<n_act; i++)
    {
      
    for(k=0;k<3;k++) 
      {
      x[ind_act[i]][k] = x_act[i][k];
      v[ind_act[i]][k] = v_act[i][k];
      } /* k */

    t[ind_act[i]] = t_act[i];
    dt[ind_act[i]] = dt_act[i];

    pot[ind_act[i]] = pot_act[i];

    for(k=0;k<3;k++) 
      {
      a[ind_act[i]][k] = a_act[i][k];
      adot[ind_act[i]][k] = adot_act[i][k];
      } /* k */
        
    } /* i */


#ifdef TIMING
  get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
  DT_ACT_CORR += (CPU_tmp_user - CPU_tmp_user0);
#endif




  /* load the new values for active particles to the GRAPE */

#ifdef TIMING
  get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

  for(j=0; j<n_act; j++) 
    {
    jjj = ind_act[j];
      
    for(k=0;k<3;k++) 
      {
      a2by18[k] = 0.0; a1by6[k] = adot_act[j][k]*over6; aby2[k] = a_act[j][k]*over2;

      xcm_tmp[k] = x_act[j][k] - xcm[k];
      vcm_tmp[k] = v_act[j][k] - vcm[k];
      } /* k */

//    g6_set_j_particle(clusterid, jjj, ind_act[j], t_act[j], dt_act[j], m_act[j], a2by18, a1by6, aby2, v_act[j], x_act[j]);
    g6_set_j_particle(clusterid, jjj, ind_act[j], t_act[j], dt_act[j], m_act[j], a2by18, a1by6, aby2, vcm_tmp, xcm_tmp);
    } /* j */

#ifdef AMD
    g6_flush_jp_buffer(clusterid);
#endif

#ifdef TIMING
  get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
  DT_ACT_LOAD += (CPU_tmp_user - CPU_tmp_user0);
#endif




#ifdef GRAPESHARING
    grape_sharing();
#endif


 
  /* Current time set to min_t */

  time_cur = min_t;

  Timesteps++;
  n_act_sum += n_act;
  n_act_distr[n_act-1]++;

/*
  for(i=1; i<N+1; i++) 
    {
    if(i == n_act) n_act_distr[i-1]++;
    }
*/


#ifdef DEBUG111
  printf("%6.0f \t %06d \t %6.0f \n", Timesteps, n_act, n_act_sum);
  fflush(stdout);
#endif




#ifdef STEVOL

#ifdef TIMING
  get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

  if(time_cur >= t_stevol)
    {
    /* Define the current mass of all star particles... */

    for(i=0; i<N; i++) m[i] = m_curr(m0[i], ZZZ[i], t0[i], time_cur);

    t_stevol += dt_stevol;
    } /* if(time_cur >= t_stevol) */
    
#ifdef TIMING
  get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
  DT_STEVOL += (CPU_tmp_user - CPU_tmp_user0);
#endif  
    
#endif




  if(time_cur >= t_disk)
    {
    diskstep++;      	

    write_snap_data();
    write_cont_data();
    
#ifdef DEBUG
    write_dbg_data();
#endif
    

#ifdef GRAPESHARINGOLD
    grape_sharing_old();
#endif 		// GRAPESHARINGOLD

    t_disk += dt_disk;
    } /* if(time_cur >= t_disk) */





  if(time_cur >= t_contr)
    {


#ifdef STARDESTR

#ifdef TIMING
  get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif

  star_destr(time_cur, &m_bh, &e_bh, &num_bh, n_act, ind_act, m, x, v, pot, pot_ext, dt);

  out = fopen("phi-GRAPE.ext","w");
  fprintf(out,"%.8E \t %.8E \n", m_bh, b_bh);
  fclose(out);

  out = fopen("mass-bh.dat","a");
  fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh);
  fclose(out);
  

    /* load the new values for particles to the GRAPE */

    nj=N; 

    /* load the nj particles to the GRAPE */

    for(j=0; j<nj; j++) 
      {
      for(k=0;k<3;k++) 
        {
        a2by18[k] = 0.0; a1by6[k] = adot[j][k]*over6; aby2[k] = a[j][k]*over2;
        } /* k */

      g6_set_j_particle(clusterid, j, ind[j], t[j], dt[j], m[j], a2by18, a1by6, aby2, v[j], x[j]);
      } /* j */

#ifdef AMD
    g6_flush_jp_buffer(clusterid);
#endif
  

#ifdef TIMING
  get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
  DT_STARDESTR += (CPU_tmp_user - CPU_tmp_user0);
#endif
    
#endif		// STARDESTR


    
    energy_contr();


#ifdef CMCORR
    for(i=0; i<N; i++)
      {
        	
      for(k=0;k<3;k++) 
        {
//        x[i][k] -= xcm[k];
//        v[i][k] -= vcm[k];
        } /* k */
        
      } /* i */
#endif



#ifdef DEF_DC
  def_dc(N, m, x, v, xcm, vcm);
#else

  for(k=0;k<3;k++)
    {
    xcm[k] = 0.0;
    vcm[k] = 0.0;
    } /* k */

    /* load the new values for particles to the GRAPE */

    nj=N; 

    /* load the nj particles to the GRAPE */

    for(j=0; j<nj; j++) 
      {
      for(k=0;k<3;k++) 
        {
        a2by18[k] = 0.0; a1by6[k] = adot[j][k]*over6; aby2[k] = a[j][k]*over2;

        xcm_tmp[k] = x[j][k] - xcm[k];
        vcm_tmp[k] = v[j][k] - vcm[k];
        } /* k */

//      g6_set_j_particle(clusterid, j, ind[j], t[j], dt[j], m[j], a2by18, a1by6, aby2, v[j], x[j]);
      g6_set_j_particle(clusterid, j, ind[j], t[j], dt[j], m[j], a2by18, a1by6, aby2, vcm_tmp, xcm_tmp);
      } /* j */

#ifdef AMD
    g6_flush_jp_buffer(clusterid);
#endif

#endif		// DEF_DC




    /* possible OUT for timing !!! */

#ifdef TIMING
    out = fopen("timing.dat","a");
    
    DT_TOT = DT_ACT_DEF  + DT_ACT_PRED + 
             DT_ACT_GRAV + DT_EXT_GRAV + 
             DT_ACT_CORR + DT_ACT_LOAD + 
             DT_STEVOL + DT_STARDISK + DT_STARDESTR;

    fprintf(out,"%.8E \t %.6E \t %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f \t %.8E %.8E %.8E \n",
            time_cur, DT_TOT, 
            100.0*DT_ACT_DEF/DT_TOT,  100.0*DT_ACT_PRED/DT_TOT, 
            100.0*DT_ACT_GRAV/DT_TOT, 100.0*DT_EXT_GRAV/DT_TOT, 
            100.0*DT_ACT_CORR/DT_TOT, 100.0*DT_ACT_LOAD/DT_TOT, 
            100.0*DT_STEVOL/DT_TOT, 100.0*DT_STARDISK/DT_TOT, 100.0*DT_STARDESTR/DT_TOT,
            CPU_time_real-CPU_time_real0, CPU_time_user-CPU_time_user0, CPU_time_syst-CPU_time_syst0);
    
    fclose(out);
#endif		// TIMING



    t_contr += dt_contr;
    } /* if(time_cur >= t_contr) */



  } /* while(time_cur <= t_end) */		// End of the main integration loop... 



  /* close the GRAPE */
  
  g6_close(clusterid);


  /* Write some output for the timestep annalize... */

  printf("\n");
  printf("Timesteps = %.0f   Total sum of integrated part. = %.0f   g6_calls on all nodes = %.0f \n", Timesteps, n_act_sum, g6_calls);
  printf("\n");
  printf("Real Speed = %.3f GFlops \t Teor Speed = %.3f GFlops \n", 57.0*N*n_act_sum/(CPU_time_user-CPU_time_user0)/1.0E+09, 57.0*N*g6_calls*48.0/(CPU_time_user-CPU_time_user0)/1.0E+09);
  fflush(stdout);


#ifdef DEBUG_extra
  tmp_file = fopen("n_act_distr.dat","w");

  printf(tmp_file,"\n");
  printf(tmp_file,"Total Timesteps = %.0f   Total sum of integrated part. = %.0f   g6_calls on all nodes = %.0f \n", Timesteps, n_act_sum, g6_calls);
  printf(tmp_file,"\n");
  printf(tmp_file,"Real Speed = %.3f GFlops \t Teor Speed = %.3f \n", 57.0*N*n_act_sum/(CPU_time_user-CPU_time_user0)/1.0E+09, 57.0*N*g6_calls*48.0/(CPU_time_user-CPU_time_user0)/1.0E+09);

  for(i=1;i<N+1;i++) 
    {
    fprintf(tmp_file,"%06d \t %.6E \t %.2f \n", i, n_act_distr[i-1], 100.0*n_act_distr[i-1]/Timesteps);
    }
  fclose(tmp_file);
#endif		// DEBUG_extra 
 

  return(0);

  }

/****************************************************************************/
/****************************************************************************/
/****************************************************************************/
